
b_train=[];
b_trainFit=[];
b1_valError=[];
b2_pred=[];
b2_corr=[];
b2_valError=[];
cellID=[];
for tnd=1:length(proj_meta(cur_site).ExpGroup)
    tnd
    Y=proj_meta(cur_site).raw_data(1,tnd).velM_raw(end-4999:end)';
    X=[];
    cnt=1;
    for gnd=1:4
        for fnd=1:length(proj_meta(cur_site).raw_data(gnd,tnd).smoothed_ROIs);
            X(:,cnt)=proj_meta(cur_site).raw_data(gnd,tnd).smoothed_ROIs(fnd).activity(end-4999:end)/median(proj_meta(cur_site).raw_data(gnd,1).smoothed_ROIs(fnd).activity);
            cnt=cnt+1;
        end
    end
    b = TreeBagger(50,X,Y,'method','r','oobpred','on','oobvarimp','on','minleaf',50);
    b_train(tnd,:)=oobPredict(b)';
    b_trainFit(tnd)=corr(b_train(tnd,:)',Y)^2;
%     b1_valError(tnd,:)=b.OOBPermutedVarDeltaError;
%     cur_vals=sort(b1_valError(tnd,:),'ascend');
%     max_vals=cur_vals(end-9:end);
%     idxvar=find(ismember(b1_valError(tnd,:),max_vals));
%     cellID(tnd,:)=idxvar;
%     b2 = TreeBagger(50,X(:,idxvar),Y,'method','r','oobpred','on','oobvarimp','on','minleaf',50);
%     b2_pred(tnd,:)=oobPredict(b2)';
%     b2_corr(tnd)=corr(b2_pred(tnd,:)',Y)^2;
%     b2_valError(tnd,:)=b2.OOBPermutedVarDeltaError;
end

figure;
for tnd=1:length(proj_meta(cur_site).ExpGroup)
    subplot(7,1,tnd)
    plot(proj_meta(cur_site).raw_data(1,tnd).velM_raw(end-4999:end)/max(proj_meta(cur_site).raw_data(1,tnd).velM_raw(end-4999:end)),'r')
    hold on
    plot(b_train(tnd,:)/max(b_train(tnd,:)),'b')
    ylim([-0.01 1])
end

figure;
for tnd=1:length(proj_meta(cur_site).ExpGroup)
    subplot(7,1,tnd)
    plot(proj_meta(cur_site).raw_data(1,tnd).velM_raw(end-4999:end)/max(proj_meta(cur_site).raw_data(1,tnd).velM_raw(end-4999:end)),'r')
    hold on
    plot(b2_pred(tnd,:)/max(b_train(tnd,:)),'b')
    ylim([-0.01 1])
end

figure;plot(b_trainFit);hold on;plot(b2_corr,'r')

all_cells=[];
for hnd=1:length(cellID)
    all_cells(hnd,:)=cellID{hnd};
end
all_cells=unique(all_cells);

figure;
set(gcf,'position',[1932         305        1898         685])
hold on
for tnd=1:length(proj_meta(cur_site).ExpGroup)
    subplot(1,length(proj_meta(cur_site).ExpGroup),tnd)
    X=[];
    cnt=1;
    for gnd=1:4
        for fnd=1:length(proj_meta(cur_site).raw_data(gnd,tnd).smoothed_ROIs);
            X(:,cnt)=proj_meta(cur_site).raw_data(gnd,tnd).smoothed_ROIs(fnd).activity(end-4999:end)/median(proj_meta(cur_site).raw_data(gnd,1).smoothed_ROIs(fnd).activity(end-4999:end));
            cnt=cnt+1;
        end
    end
    imagesc(X(:,all_cells)')
    set(gca,'clim',[1 5.5])
end
figure;
set(gcf,'position',[1930          26        1897         216])
hold on
for tnd=1:length(proj_meta(cur_site).ExpGroup)
    subplot(1,length(proj_meta(cur_site).ExpGroup),tnd)
    plot(proj_meta(cur_site).raw_data(1,tnd).velM_raw(end-4999:end)/max(proj_meta(cur_site).raw_data(1,tnd).velM_raw(end-4999:end)))
    ylim([0 1])
end

figure;
set(gcf,'position',[1932         305        1898         685])
for tnd=1:length(proj_meta(cur_site).ExpGroup)
    subplot(1,length(proj_meta(cur_site).ExpGroup),tnd)
    X=[];
    cnt=1;
    for gnd=1:4
        for fnd=1:length(proj_meta(cur_site).raw_data(gnd,tnd).smoothed_ROIs);
            X(:,cnt)=proj_meta(cur_site).raw_data(gnd,tnd).smoothed_ROIs(fnd).activity(end-4999:end)/median(proj_meta(cur_site).raw_data(gnd,1).smoothed_ROIs(fnd).activity(end-4999:end));
            cnt=cnt+1;
        end
    end
    for knd=1:length(all_cells)
        hold on
        plot(X(:,all_cells(knd))'-3*(knd-1))
    end
%     set(gca,'clim',[1 5.5])
end
figure;
set(gcf,'position',[1930          26        1897         216])
hold on
for tnd=1:length(proj_meta(cur_site).ExpGroup)
    subplot(1,length(proj_meta(cur_site).ExpGroup),tnd)
    plot(proj_meta(cur_site).raw_data(1,tnd).velM_raw(end-4999:end)/max(proj_meta(cur_site).raw_data(1,tnd).velM_raw(end-4999:end)))
    ylim([0 1])
end








win_size=900;
% find playback onset cells
all_wind=[];
for tnd=1:length(proj_meta(cur_site).ExpGroup)
    if length(proj_meta(cur_site).raw_data(1,tnd).nbr_frames)>5
        order=[];
        for xnd=1:length(proj_meta(cur_site).raw_data(1,tnd).nbr_frames)
            if proj_meta(cur_site).raw_data(1,tnd).nbr_frames(xnd)==5000
                if xnd==length(proj_meta(cur_site).raw_data(1,tnd).nbr_frames)
                    order(xnd)=3;
                else
                    order(xnd)=1;
                end
            else
                order(xnd)=2;
            end
        end
    else
        order=[1 2 1 2 3];
    end
    a=1:length(order);
    order(max(a(order==2))+1:end)=3;
    exp_frames=ones(2,length(order));
    for bnd=1:length(order)
        if bnd>1
            exp_frames(1,bnd)=exp_frames(2,bnd-1)+1;
        end
        exp_frames(2,bnd)=sum(proj_meta(cur_site).raw_data(1,tnd).nbr_frames(1:bnd));
    end
    play_ind=find(order==2);
    wind_ind=[];
    for vnd=1:length(play_ind)
        for pnd=exp_frames(1,play_ind(vnd)):exp_frames(2,play_ind(vnd))-win_size
            if sum(proj_meta(cur_site).raw_data(1,tnd).velM_ind(pnd:pnd+win_size))==0 && sum(proj_meta(cur_site).raw_data(1,tnd).velP_ind(pnd:pnd+win_size))>win_size/3
                wind_ind=pnd;
                break
            end
        end
    end
    if isempty(wind_ind)
        error('Decrease window size')
    else
        all_wind(tnd)=wind_ind;
    end
end

figure;
set(gcf,'position',[1930          26        1897         216])
hold on
for tnd=1:length(proj_meta(cur_site).ExpGroup)
    subplot(1,length(proj_meta(cur_site).ExpGroup),tnd)
    plot(proj_meta(cur_site).raw_data(1,tnd).velP_raw(all_wind(tnd):all_wind(tnd)+(win_size-1)),'r')
    filt=ftfil(proj_meta(cur_site).raw_data(1,tnd).velP_raw(all_wind(tnd):all_wind(tnd)+(win_size-1)),30,0,10);
    hold on
    plot(smooth2(filt,10),'b')
    plot(proj_meta(cur_site).raw_data(1,tnd).velM_raw(all_wind(tnd):all_wind(tnd)+(win_size-1)),'g')
    xlim([0 win_size])
end

% cells=[15 66 17 69 82 99 104 155 215 226 245 251 334 335 388 419];
cells=[2 12 13 15 17 37 52 66 80 85 86 155 168 203 212 223 224 226 230 264 303 317 401 419 427 429 444 477 482 562];

% b_train=[];
% b_trainFit=[];
% trainP=[];
% b_pred=[];
% b_predFit=[];
% predP=[];

b_train=[];
b_trainFit=[];
trainP=[];
b_pred={};
b_predFit=[];
predP={};
% for lnd=1:5
for tnd=1:length(proj_meta(cur_site).ExpGroup)
%     if length(proj_meta(cur_site).raw_data(1,tnd).nbr_frames)>5
%         order=[];
%         for xnd=1:length(proj_meta(cur_site).raw_data(1,tnd).nbr_frames)
%             if proj_meta(cur_site).raw_data(1,tnd).nbr_frames(xnd)==5000
%                 if xnd==length(proj_meta(cur_site).raw_data(1,tnd).nbr_frames)
%                     order(xnd)=3;
%                 else
%                     order(xnd)=1;
%                 end
%             else
%                 order(xnd)=2;
%             end
%         end
%     else
%         order=[1 2 1 2 3];
%     end
%     a=1:length(order);
%     order(max(a(order==2))+1:end)=3;
%     exp_frames=ones(2,length(order));
%     for bnd=1:length(order)
%         if bnd>1
%             exp_frames(1,bnd)=exp_frames(2,bnd-1)+1;
%         end
%         exp_frames(2,bnd)=sum(proj_meta(cur_site).raw_data(1,tnd).nbr_frames(1:bnd));
%     end
    tnd
    filt=ftfil(proj_meta(cur_site).raw_data(1,tnd).velP_raw,30,0,10);
    runP=smooth2(filt,3);
    %     Y=runP(all_wind(tnd)-4:all_wind(tnd)+(win_size-1)-4)';
    Y=runP(cor_snips{tnd})';
    %     Y=runP([fb_frames(1,1):fb_frames(2,1) fb_frames(1,2):fb_frames(2,2)])';
    data=[];
    cnt=1;
    for gnd=1:4
        for fnd=1:length(proj_meta(cur_site).raw_data(gnd,tnd).smoothed_ROIs);
            raw=proj_meta(cur_site).raw_data(gnd,tnd).smoothed_ROIs(fnd).activity;
            data(:,cnt)=raw/median(raw);
            cnt=cnt+1;
        end
    end
    %     X=data(all_wind(tnd):all_wind(tnd)+(win_size-1),cells);
    X=data(cor_snips{tnd},cells);
    %     X=data([fb_frames(1,1):fb_frames(2,1) fb_frames(1,2):fb_frames(2,2)],cells);
    %     X=X(5-(lnd-1):end,:);
    %     Y=Y(1:end-(5-lnd));
    b = TreeBagger(50,X,Y,'method','r','oobpred','on','oobvarimp','on','minleaf',5);
    b_train(tnd,:)=predict(b,X);
    trainP(tnd,:)=Y;
    b_trainFit(tnd)=corr(b_train(tnd,:)',trainP(tnd,:)')^2;
    %
    %         fb_frames=exp_frames(:,order==1);
    %         b_pred(tnd,:)=predict(b,data([fb_frames(1,1):fb_frames(2,1) fb_frames(1,2):fb_frames(2,2)],cells));
    %         predP(tnd,:)=runP([fb_frames(1,1):fb_frames(2,1) fb_frames(1,2):fb_frames(2,2)]);
    %         b_predFit(tnd,:)=corr(b_pred(tnd,:)',predP(tnd,:)')^2;
    %     trainP{tnd}=Y;
    %     for knd=1:size(data,2)
    %         b_train{tnd}(knd,:)=predict(b,X(:,knd));
    %         b_trainFit(knd,tnd)=corr(b_train{tnd}(knd,:),trainP{tnd})^2;
    %     end
    %     fb_frames=exp_frames(:,order==1);
    %     b_pred{tnd}=predict(b,data([fb_frames(1,1):fb_frames(2,1) fb_frames(1,2):fb_frames(2,2)],cells));
    %     predP{tnd}=runP([fb_frames(1,1):fb_frames(2,1) fb_frames(1,2):fb_frames(2,2)]);
    %     b_predFit(tnd,:)=corr(b_pred{tnd},predP{tnd}')^2;
end
% end

figure;
for tnd=1:length(proj_meta(cur_site).ExpGroup)
    subplot(length(proj_meta(cur_site).ExpGroup),1,tnd)
    plot(trainP(tnd,:),'r')
    hold on
    plot(b_train(tnd,:),'b')
end
figure;plot(b_trainFit);

figure;
for tnd=1:length(proj_meta(cur_site).ExpGroup)
    subplot(length(proj_meta(cur_site).ExpGroup),1,tnd)
    plot(predP{tnd},'r')
    hold on
    plot(b_pred{tnd},'b')
end
figure;plot(b_predFit);


figure;
for tnd=1:length(proj_meta(cur_site).ExpGroup)
    subplot(length(proj_meta(cur_site).ExpGroup),1,tnd)
    plot(predP(tnd,:),'r')
    hold on
    plot(b_pred(tnd,:),'b')
end

figure;plot(b_predFit);